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ABSTRACT 



Context. Recent results from interferometry and asteroseismology require models of rapidly rotating stars that are more and more 
precise. 

Aims. We describe the basic structure and the hydrodynamics of a fully radiative star as a preliminary step towards more realistic 
models of rotating stars. 

Methods. We consider a solar mass of perfect gas enclosed in a spherical container. The gas is self-gravitating and rotating, and is the 
seat of nuclear heating, and heat diffusion is due to radiative diffusion with Kramers type opacities. Equations are solved numerically 
with spectral methods in two dimensions with a radial Gauss-Lobatto grid and spherical harmonics. 

Results. We computed the centrifugally flattened structure of such a star: the von Zeipel model, which says that the energy flux is 
proportional to the local effective gravity is tested. We show that it overestimates the ratio of the polar to the equatorial energy flux by 
almost a factor 2. We also determine the Brunt-Vaisala frequency distribution and show that outer equatorial regions in a radiative zone 
are convectively unstable when the rotation is fast enough. We compute the differential rotation and meridional circulation stemming 
from the baroclinicity of the star and show that, in such radiative zones, equatorial regions rotate faster than polar ones. The surface 
differential rotation is also shown to reach a universal profile when rotation is slow enough (less than 36% of the breakup one), as 
long as viscosity and Prandlt number remain small. 
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1. Introduction 

Recent observations of nearby stars with optical or infrared inter- 
ferometers have opened a new window on rapidly rotating stars. 
The w ork of Domiciano de Souza et al. (2005), Zorec et 
(l2005h. McAlister et al. (2005), Peterson et al. (2006^ 
iPeterson et al. (2006b) . K erv ella & Domician o de Souza , ( 2006h 
and lAufdenberg eralj ( 20061) show that observational tech- 
niques are now able to constrain the surface distribution of 
the energy flux of these stars along with their precise shape 
on the background sky. Data are usually fit by a simple 
brightness law derived from a R oche-von Zeipel model (e.g. 
iDomiciano de Souza et all l2002h where it is assumed that the 
shape of the star is determined by a Roche model (all the mass 
is concentrated in the center), while the brightness distribution 
obeys the von Zeipel law, namely the energy flux is proportional 
to local effective gravity. If such a model is efficient at giving 
a first interpretation of the observation and especially getting 
an idea of the orientation of the stellar rotation axis, as well as 
its angular velocity, one now wishes to go beyond this simple 
model. This step is difficult, however, as it means building a 
self-consistent two-dimensional model of a star and then being 
able to make it evolve till the assumed age of the observed 
object. 

Many attempts have been made to reach this goal, but most 
of them have not gone beyond the construction of barotropic 
stars without any dynamics or evolution (see the recent work 
oflRoxburgh, 2004, Jackson et al., 2005). The inclusion of fluid 
dynamics is, however, important for the models to be physi- 
cally self-consistent. Rotating stars are indeed baroclinic, and 



fluid flows are present everywhere, even in radiative zones. Their 
long-term effect is the r nixing of ele ments and the transport of 
angular momentum (see lZahiull992l) . so they cannot be ignored 
in the evolution of rotating stars. 

Previous work on models of baroclinic stars is rare. Only 
lUrvu & Eriguchil d 1 994 1 1 9951) constructed global 2D models of 
rapidly rotating stars with baroclinicity, but their neglect of vis- 
cosity removed any meridional circula tion and left d ifferential 
rotation undetermined. More recently, Rieutord] ( 120061) used a 
Boussinesq model (i.e. a nearly incompressible fluid) for investi- 
gating the baroclinic flows on a global scale. In this case, it could 
be shown how viscosity determines the differential rotation and 
that this rotation strongly depends on the Prandtl number. 

The present work aims at describing a more realistic, physi- 
cally self-consistent, approach of a rotating, fully radiative star. 
Of course this is a model and real physics is still very simpli- 
fied, but no arbitrary assumption needs to be made. In a few 
words, our model can be described as follows: we consider a 
perfect gas contained in a rigid sphere; the gas is rotating and 
self-gravitating; the mass inside the sphere is high enough so that 
the central temperature and density permit nuclear reactions. The 
gas is a viscous fluid conducting heat through radiative diffusiv- 
ity derived from Kramers type opacities. The bounding sphere 
allows us to impose simple (but physically relevant) boundary 
conditions on the gravitational potential, temperature, pressure, 
and velocity field. Roughly speaking, these conditions specify 
(i) that the gravitational potential connects to the one pervading 
vacuum, (ii) that temperature connects to the temperature field 
of a medium with a constant absorption coefficient radiating like 
a black body, (iii) that some radial stress can be accepted by 
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the bounding sphere, but (iv) that matter slips freely and tangen- 
tially on the sphere. This set-up, which is similar to a laboratory 
experiment, simplifies the numerical resolution of the equations 
by allowing the use of spherical coordinates while the distribu- 
tion of mass is not spherically symmetric. Fortunately, the bulk 
properties of this "star" do not depend very much on the physics 
imposed at the outer boundary. The next step will be to replace 
the sphere by a spheroid that fits an isobar perfectly; however, 
this step, which needs to use spheroidal coordinates with multi- 
domains, is numerically challenging and can be undertaken only 
on a physically understood system. 

We have voluntarily removed the convection zones that need 
a generalization of the mixing-length theory in two dimensions 
including the effects of rotation. Here we focus on a rotating 
radiative star and examine the dynamical properties of its baro- 
clinic state. 

The paper is organized as follows. We next describe our 
model (Sect. 2) and the numerical methods we use for solv- 
ing the equations (Sect. 3). We then describe some tests of the 
model showing internal accuracy and how they compare to one- 
dimensional models (Sect. 4). In Sect. 5 results are presented and 
discussed, and conclusions follow. 



2. The model 

We consider a viscous, self-gravitating fluid rotating with a mean 
angular velocity Q - fiC;. The fluid is supposed to be inside 
an immutable sphere of radius R, so the outer boundary can be 
taken at r = This simplifies the problem by allowing spheri- 
cal coordinates to represent the variables, although the external 
boundary does not fit an isobar 

In a steady configuration, the fluid is governed by the equa- 
tions 



= AnGp 
pTv ■ Vi = -V ■ f + e 



(1) 



p(2nxv-Hv- Vv) = -VP-pV(0- -Q'-r'- sin^0)-t-f„ 



V ■ (pv) = 



where ^ is the gravitational potential, p the density, s the specific 
entropy, T the temperature, v the fluid velocity with respect to 
a frame rotating at angular velocity £1. Also, P is the pressure 
and s the nuclear energy generation rate per unit volume. For 
a compressible fluid with constant dynamical viscosity yU, the 
viscous force can be written as 



(2) 



For the energy flux, we will consider only radiative transport, 

F = -xVT (3) 

wheTexiP' T) is the radiative difFusivity. The equation of state is 
that of an ideal gas mixed with photons; hence, we use 



(4) 



where I^m is the gas constant divided by the mean molecular 
weight and a the radiation density constant. 



2. 1 . Stellar microphysics 

We assume that the fluid is in the conditions of a stellar plasma 
typical of the sun's radiative zone. The nuclear energy generation 
rate per unit volume has the form 

s = soXVt-'"c-'^'"' (5) 

where X is the hydrogen mass fraction. Fo r the constants eo and b 
we adopt the values of the CESAM code dMorell [19971) . namely 

b = 3600 and so = 8.37 lO'" (cgs) (6) 

for the pp-chain. 

In order to describe the radiative transport of energy, we use 
the opacity given by power laws: 



so that the radiative difFusivity can be written as 
16o-T^ 



X = 



3kp 



(7) 



(8) 



In particular we use [3 = 1.97541, ri = 0.138316, and 
kq = 7.1548412 10'^ cgs, following some fitting formulae^] pro- 
posed by J. Christensen-Dalsgaard for simple solar models (e.g. 
IChristensen-Dalsgaard & Reiteil[T995l) . 



2.2. Boundary conditions 

In order to solve the full system of equations, we must complete 
it with boundary conditions. At the center of the star (r - 0) 
we just need to impose regularity on the solutions. At the outer 
boundary, i.e. on the rigid bounding sphere, we impose that the 
fluid slips freely and thus use stress-free boundary conditions. 
Hence, we have 



Ur - (Ti-e - (Trtp = at r - R 



(9) 



where cr,.g and cr^^ are the horizontal components of the stress 
on the outer sphere. 

Since the sphere is rigid, it can support some normal stress, 
and pressure is not constant on it any more than the density and 
temperature. As the true surface of the star is outside the box, 
pressure needs to be fixed somewhere. We thus impose the polar 
pressure. This condition imposes the place where our container 
"cuts" the star envelope. 

The usual boundary conditions on the temperature are that 
the boundary radiates like a black body. We thus write 

dT 4 
or 

Using the expression of radiative diffusivity, this condition is 
equivalent to 

dT 

— +(TtT^O 
or 

where o-j - 3pK/16. Obviously, o-j depends on latitude, but to 
avoid unnecessary complications (nonlinearities) we take it as 
a constant. Physically, this means that at the boundary the fluid 
meets a medium with constant absorption coefficient radiating 
like a black body. 

Finally, the boundary conditions on the gravitational poten- 
tial are that it matches the vacuum solution, which vanishes at 
infinity. 



' The proposed formula is ^ = 7 + 7- with /c, as in Q and 
= 1.6236784 io-33p0.407895 7-9.28289. J^iy jfjg "interior" opacity, 

leaving the "envelope" one. 
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2.3. Dimensionless equations 

We scale the physical quantities, so that the solutions are 
of order unity. For that, we choose the free-fall time scale 
A - (47rGpc) the radius of the container/? as the length scale, 
the central values of temperature and density pc for the corre- 
sponding scales, and the gas constant "Rm for the entropy scale. 
With dimensionless variables, the equations are now written as 

r 

■^pTu ■ Vi = V ■ (x^T) + As 

p (2Q X M + M ■ Vm) = -Vp - pV^eff + EFu (10) 

V ■ (pM) = 

[p^nJ{pT+^T') 

where 0eff - <p- j^'^r^ sin^ 6 and F„ - Au + | V(V ■ «). We also 
introduced the dimensionless numbers 



E = 



pA 



Xc 



Here, £ is a kin d of Ekm an number associated with the limit an- 
gular velocity ^AnGpc, and P the Prandtl number at the center. 
The Ekman number measures the importance of viscous force 
relative to Coriolis force, and its true definition is Ejil; here 
we had to keep out the non-dimensional rotation rate for con- 
venience. In stellar conditions EjQ. <s: 1. Other, more stellar, 
dimensionless numbers have been introduced: 



'RmJcA 



A 



XcTc 



and Pc 



Kr + 



aT*A- 



where tTc is the dimensionless central gas pressure and /3c the 
usual ratio of gas to total pressure at the star center Once we fix / ~ / , fM^i (^) 
the microphysics of the problem, the dimensionless numbers tTc 
and A should appear as part of the solution, as they determine the 
central temperature and density (or mass) of the configuration. 

The quantities x and s have been scaled using their central 
values Xc and Sc, so that the dimensionless radiative diffusivity 
and energy generation rate are given by 



3. Numerical method 

One-dimensional stellar evolution codes usually use a few thou- 
sand of grid points to model and follow the evolution of a star 
With the same kind of discretization (i.e. finite differences), the 
number of grid points would grow by a factor of a few hundred 
for a two-dimensional model. This large number of points is cer- 
tainly very heavy to manage; two-dimensional models require 
more efficient discretization, which is why we adopted spectral 
methods. 

Spectral methods are best-known for t heir pr e cision 
when the number of grid points i s hmi ted (iPevrell l2002l 
iGrandclemeiiH l2006i iBonazzola etaLl Il999l) . On the dark side, 
they have never been used in stellar modeling (to our knowledge) 
and are rather "touchy" as far as stability is concerned; more par- 
ticularly, discontinuities require special care like a multi-domain 
approach does. Nevertheless, the precision they can achieve is a 
major advantage, especially when the oscillation spectrum of a 
star needs to be computed. The advances in observational aster- 
oseismology, with s pace missions like COROT, indeed require 
very precise models jReese et al.Ll2006h . 

The discretization of the equations is made via an expansion 
of the unknowns onto the spherical harmonic basis for the angu- 
lar part and using Gauss-Lobatto collocation nodes for the radial 
part. This latter grid is associated with C hebyshev polynomials 
(e.g. lCanuto et aUl2006irFornbergill998l) . 

3. 1 . Projection on the sphericai harmonics 

The variables are thus first projected on the spherical harmonics. 
Hence, for the scalar quantities: 



(15) 



Note that, due to the symmetry of the problem, solutions are 
axisymmetric and therefore m - 0. They are also symmetric 
with respect to the equatorial plane, so scalar quantities need 
only s pherical harmon ics of an even order. For the velocity, we 
follow lRieutordl (Il987l) and use the expansion: 



X ■ 



^Tf^^3p-v-i and £ = p2y,-2/3g-„(r->'^-i) 



P 



— 1 /"^ 

where a - bTc . Using these expressions and writing the en- 
tropy gradient as a function of the temperature and pressure gra- 
dients, the equation of energy balance reads 



(11) H = ^ Eui{r)R^, + Ev,(r)S'^ + wi(r)T*l, 



(16) 



where the vectorial spherical harmonics are 



R 



(12) 



(13) 



r p''+^ 

Ar H- V In y ■ Vr a ■ M = -A „ 

where we introduced the vector 

^-^-^[Tz-r^'Vc^p^ 

which makes explicit the entropy gradient; y is the adiabatic ex- 
ponent. 

The scal i ng of t he velocity field requires some consideration. 
In iRieutordI (120061) it is shown that the meridional circulation, 
i.e. the radial and meridional components of the velocity field, is 
proportional to the Ekman number; therefore, we write: 



dY"' 



00 sm6 dip 



1 dY' 



dY, 



(17) 



' sine 5^ ' do ''- 

The projection of the momentum equation yields three equa- 
tions: 



^+2<p)Q[(/-l)a,w,_i 
dr 



-E' 



A-dFui 8 dui 
3 dr^ 3r dr 



- (^ + 2)a/+iw,+i] 
(/(/ -I- 1) H- 8/3) 

Ul 



U - EUrCr + EugCe + "inC^ 



(14) 



/(/+l)dv, 

+ — — — Vl 



3r dr 



3r2 



= rhsu, 



(18) 



The system (fTOl l is the set of nonlinear PDE that governs the 
radiative zone of a rotating star when no turbulence or magnetic 
field troubles the place. It is this system that will give us the first 
view of the dynamics inside the central parts of a rapidly rotating 
fully radiative star. 



— +2{p)D. 
r 



l-l 1 + 2 

— — a/vv;-! -I- - — -ai+iwi+i 
I l+l 



d\ 2 dv, 41(1 + 1) 
dr^ r dr 3r^ ' 



1 dut 8 
3r dr 3r^ ' 



rhsv; (19) 
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3.3. Algorithm for iterations 



2<p>Q 



— ui-i - - — - 
I l+l 



l-l 1+2 

aivi-i - - — -ai+ii'i+i 



d^wi 2 dwi 1(1+1) 



-Wl 



rhsw/, (20) 



dr^ r dr 

where we have introduced the spherically averaged quantities 
Mr) 



and the coupling coefficient 
/ 



(21) 



V4/2 - 1 
The other equations yield 



(22) 



<f>i ^Pi 



d^4>i 2d4>i 1(1+1) 
dr^ r dr 

dui 2m/ 1(1 + 1) 

— + 5 — V/ = rhsd, 

dr r 



d^T,l2 d \dT, 1(1+1) 



(23) 



(24) 



Ti-r (fl,) ui = rhst, (25) 



The right-hand terms are essentially non-linear terms that are 
gathered in Table 1. In their expressions 6 denotes the non- 
spherical part of a variable, i.e. 5 f - f - (/). 

3.2. Boundary conditions 

Boundary conditions have to be expressed on the radial func- 
tions. Regularity of the solutions at the star center demands the 
regularity of the radial functions at r = 0. As far as the pressure is 
concerned, we fix the polar value to which fixes the constant 
needed by the solution of the momentum equation . Stress-fre e 
boundary conditions on the velocity field give (iRieutordill987h 



(ui = 

dvi 

r- y; = 

dr 

dwi 

r— w, = 

dr 



(29) 



at the surface. The gravitational potential at the surface matches 
that of the vacuum 

r— H-(/-H 1)0, = at r=l. (30) 
dr 

Finally, the radial functions of the temperature verify 



r'El+o-TTi^O (31) 
dr 



Equations are solved iteratively because of n onlinearities. We 
follow and generalize the algorithm used by iBonazzola et al.l 
(1998) for polytropic stellar models. The idea is to set the equa- 
tions in the form 

■Cn(yn+\) = RHS(y„,y„-i) 

where X.n is a linear operator that may depend on the n-th iter- 
ate. Such a scheme is also known as the fixed-point algorithm. 
Convergence depends on the norm of the linear operator and the 
closeness to the solution; the dependence of RHS ony„-i stands 
for the relaxation terms. Unlike Henyey's method, which is a 
Newton type algorithm, this scheme does not need the computa- 
tion of the Jacobian matrix, which is quite convenient with our 
discretization and high order operators. 



4. Tests of the model 

4.1. Internal accuracy 

In Fig. [T] we show the spectra of the different physical quanti- 
ties. As these spectra are scaled by the maximum value, they 
show the internal relative precision of the solutions due to trun- 
cation. In this example, spectral convergence is very good on the 
gravitational potential, pressure, temperature and less good on 
the velocity components; this last quantity is indeed facing rapid 
variations in the boundary (Ekman) layers. 

Although spectra give a good idea of the precision of the so- 
lution, they do not tell the whole story. Indeed, round-off eiTors 
also affect the solution as they are usually amplified by the reso- 
lution of linear s ystems (backward error) and may have a devas- 
tating effect (e.g. lValdettaro et an,l2007l) . In our case, round-off 
errors have been tested on "eigenvalues" like and A and were 
found to be below 5 10 



4.2. The virial equality 

A standard test of the accuracy of the computations of 
a rotating, self-gravita ting flui d is the virial e quality (see 
lOstriker & Bodenheime r, 1968, Eriguchi & Miilleii ri985). This 
is indeed an integral form of the momentum equation, along with 
mass conservation. Since this equality has only been used with 
barotropic stars, we give its full expression here as is appropriate 
for a baroclinic star 

We first rewrite the momentum and continuity equations as 



(2£l A pu)i + pujdjUi - -pdi(p + pD}s(es)i + djaij (32) 
djpuj = 0, (33) 

where cr,y stands for the stress tensor (including the pressure), 
namely: 



a-ij = -P6ij + n \diVj + djVi - -(dkVk)5ij^ . 

On the outer boundary, we imposed stress-free boundary 
conditions, i.e. - o-,g = cr,.^ - 0. The volume integral of 
the scalar product of ( [32] ) with r yields the virial identity 



at the surface. 



2rrei + 10} + W + 3P + I, + 2Q.L, = 0, 
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RHS„,o„, = ^ (rhsu,/f? + rhsv,S° + rhsw/T^) 



rhst : 



a 

or 

Pi. 

rde 



-26pQ. sin + uy cos ff) - p \ — {u,- + ug cot 9) + u, — — t- 

r or 



"'■57- 


Ug dUy \ 


T'dej 


/ KflK,- dllg 
+ U,— + 

\ r or 


llg \ 









er+ 



eg+ 



rhsd = -M • V Inp 

/;+3 



^^^„,3."(-"M + |5«.«-V51n^.Vr 



(26) 

(27) 
(28) 



Table 1. Non-linear terms that appear in the right-hand sides of the equations of motion ( I18H20I ). ( l24b and ( |25] ). 




^5 




Fig. 1. Spectral convergence of the various physical quantities of the model. On left: scaled coefficient of the Chebyshev polyno- 
mial expansion (maximized over the spherical harmonic degrees); on right: the corresponding spectral coefficients of the spherical 
harmonics maximized over the Chebyshev coefficients. These curves have been obtained with a model where ps - 10"^, E = 10"^, 
Q = 0.07, a = 15.2, o-T = 1.5. 



where we introduced 



as the surface stress and 



2X; 



Tiei - — I pu^dV 

hv) 



as the kinetic energy in the rotating frame. 



1 pr^sin^OdV 

J(V) 

as the momentum of inertia, 
W ■ 



5/' 

^ J(V) 



pcpd^r 



as the gravitational potential energy. 



• I cr,/c 

J(V) 



3P = - o-iidV 

J(V) 

as the internal energy, 
L = I no-ijdS j 



Jr A pudV 
(V) 



as the relative angular momentum z-component. The last two 
integrals may be vanishingly small if the background rotation 
contains all the angular momentum and if the boundary of the 
domain is near the zero pressure surface. 

We give in Table |2] the value of the virial equality scaled, 
as usual, by the gravitational energy W. We have used 32 radial 
points, Lmax = 16, Ekman number E - 10"^ and surface pres- 
sure Ps - 10"^. For comparison, the break-up velocity for these 
models is Q.ic ~ 0.2. The virial equality is ver ified with errors 
less th an 10"^ thanks to spectral precision as in'Bonazzola et alj 
Jl998|). These val i ies can be co mpared with the lO '"* obtained by 
lUrvu & Eriguchj] (11994 119951) or 10 ^ - 10"^ bv lJackson etatl 
(I2005i) resulting from the finite difference discretizations. 
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n 


Virial 


lo-** 


-5.4 10^** 


0.05 


-1.4 10^** 


0.1 


-1.310-** 


0.15 


1.4 10-' 



Table 2. Values of the virial equality scaled by the gravitational 
energy for different rotation velocities. 





ESTER 


CESAM 


Mass 


3.05 


3.00 


Radius 


2.40 


2.29 


Luminosity 


107. 


107. 


Pc 


129 


132 




1.3910-3 


1.48 10-3 


A 


232 


219 



Table 3. Comparison of the results between the two-dimensional 
code ESTER run at zero rotation and the one-dimensional code 
CESAM. 



4.3. Calibration w/fh one-dimensional models 

In the dimensionless formulation, the two numbers A and tTc are 
of special importance. They give the radius, the central temper- 
ature, and the central density once the mass is fixed. tTc is indeed 
a non-dimensional measure of the central pressure while A of 
central heating. Using the expression of these numbers we note 
that their product 

'Re, 

AtTc = 

47tGPcXc 

depends only on central density and central temperature; thus, its 
determination yields a first relation between these two quantities. 
The expression of the mass 

M^PcR^ I pdV 

J(V) 

combined with that of A gives another relation between Tc and 
Pc, thereby allowing for the determination of R, T^Pc- Actually, 
for convenience we prefer fixing the central temperature instead 
of the mass and deriving the other quantities. 

As a first test of the model, we compared the results ob- 
tained for a non-rotating configur ation with th e results of a one- 
dimensional code like CESAM (Morell lT997h . We show the re- 
sults in Table[3] There, we used a central temperature of 3.05 10^ 
and a hydrogen mass fraction of X=0.712. The physics of the 
one-dimensional code was made as close to ours as possible, 
but surface boundary conditions, implying some superficial con- 
vection, still make noticeable differences. (The three-solar mass 
was chosen to minimize the effects of convection although the 
physics does not strictly apply in this case.) Nevertheless, the 
results displayed in this table show that our model is close to 
models computed with traditional stellar evolution codes. Thus, 
we have a self-gravitating fluid in our "box", close to the stellar 
conditions. 

We do not ask for more precision at the moment since all 
the physics is not implemented; we leave for future work a more 
detailed comparison. Presently, we estimate that mass and tem- 
perature distributions are close enough to "reality" that our fluid 
flows are meaningful (as long as turbulence is negligible). 




Fig. 2. Isobars for a model with Q. - 0.07. An equatorial 
Keplerian velocity for this case is reached when = 0.082, 
so the star is rotating at 85% of its break-up velocity. Other pa- 
rameters are standard (see text). The thick isobar underlines the 
last complete isobar allowed by the container. 



0.30 



0.25 - 



0.80 - 



•V 0.15 - 



0.10 - 




0.05 - 



0.00 I , - -I--,-" , , \ , , , \ , , , \ , , , 

0.0 0.2 0.4 0.6 0.8 1.0 

r 

Fig. 3. Flatness e = 1 - tt^ of the isobars as a function of the 

equatorial radius for our model (solid line) and for the polytrope 
of index 4.37 rotating at 82% of the critical velocity, ps - 10^^, 
E = 10-^ Q. = 0.07, a = 15.2, o-j = 1.5. 



5. Results 

For all the examples shown below, we have used a polar pres- 
sure of ps = 10-^ (scaled by central pressure), a rotation rate 
of Q = 0.07, which is about 82% of the critical angular veloc- 
ity, a central temperature of Tc = {b/a = 15.2)-' - 1.33 10^ K, 
cr J- = 1.5 andySc = 1 (we neglect radiation pressure). The Ekman 
number was set to 10-^ and Prandtl number to zero. The mass 
enclosed in the container was close to one solar mass. 
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1.2 




0.2- 
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Fig. 4. Normal radiative flux at the last isobar (p = /?eq), real 
values (solid line) and the von Zeipel model (dotted line). The 
von Zeipel model is calculated from the normal effective gravity; 
both curves have been scaled by their polar value. 



Table m summarizes the properties of the models when rota- 
tion is increased. Numerics use a radial grid with Nr=64 points 
and the maximum degree of spherical harmonics is Ln,ax=32. 
We use p, = lQ-\ Tc = 1.3285 10^ K, or = 1.5 and yS, = 1; 
hydrogen mass fraction is X=0.71 and Z-Q. Ekman number is 
E = 10"^. We note that the mass increases slightly, as well as 
the luminosity, while the radius and the central density decrease 
slightly. This is a consequence of our choice to keep the central 
temperature and the ratio of polar pressure to central pressure 
constant. 

5.1. Internal structure 




Fig. 6. Square of the Brunt- Vaisala frequency A^^ when the rota- 
tion is n = 0.03. 



The first view of the internal structure is given by the pressure 
distribution as shown in Fig. |2] We emphasize the "last isobar" 
with a thick line. This is the last complete isobar fitting in the 
container, and we use it below to appreciate the latitude varia- 
tions of the flux and temperature. 

Mass distribution can be appreciated by the flatness profile 
of the isobars, as shown in Fig.|3] It is well known that a stellar 
envelope with power law opacities behaves like a polytropic en- 
velope of index n - {J3 + 3)/{ri + 1), which is ~ 4.37 in our case. 
Thus, for comparison, we also show the flatness of the isobar of 
a fully polytropic star rotating at 85% of the breakup velocity. 
The curves show that the polytrope is much more centrally con- 
densed than the star, as the flatness of isobars are almost zero at 
the center 

We then computed the energy flux surface density at the 
last isobar contained in the domain (Fig. |2]) and compared it 
to the one given by the von Zeipel model. We recall that this 
model assumes that the star is barotropic, so that density and 
temperature depend only on the effective potential {i.e. gravita- 
tional plus centrifugal); it follows that the energy flux is propor- 
tional to the local effective gravity. Figure |4] shows that the von 
Zeipel model overestimates the contrast of the polar and equato- 
rial flux. Actually, the ratio of polar to equatorial flux is 3.68 for 
our model, while it is 6.2 for von Zeipel one, almost two times 
higher. 



As far as temperature is concerned, we computed the varia- 
tion of temperature along the last isobar As shown in Fig.|5] the 
temperature drops noticeably at the equator, i.e. by a few 10^ K. 
The temperature drop shows that isothermal surfaces are more 
spherical than the isobars. 

Next w e considered the local Brunt-Vaisala frequencywhich, 
as shown in lRieutordI (l2006h . controls the baroclinic torque. This 
quantity is plotted in Fig. |6]for a rotation of Q = 0.03 and in 
Fig. [Tjfor the usual rotation D. - 0.07. Quite interestingly, we 
see that if rotation is fast enough, the local Brunt-Vaisala fre- 
quency is imaginary in the equatorial region revealing the ap- 
pearance of a convection zone driven by the cool equator. Of 
course in such a region the model that we have used is no longer 
valid as one should take into account the entropy mixing of tur- 
bulent convection. However, we think that we show here, for 
the first time, the emergence of an equatorial convection zone 
due a fast rotation. In Fig. [8] we underline the anisotropy of the 
squared Brunt-Vaisala frequency distribution by comparing the 
poleward and equatorward profiles. This anisotropy will impact 
on the low-frequency spectrum of oscillations (gravito-iner tial 
modes) as can be surmised from lDintrans & Rieutordl(l2000h . 
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Mass riVT^^ 


1.003 


1.009 


1.018 


1.032 


1.050 


1.074 


1.105 


Radius (Rq) 


0.9955 


0.9941 


0.9918 


0.9883 


0.9833 


0.9765 


0.9674 


Luminosity (Lq) 


0.803 


0.805 


0.808 


0.811 


0.816 


0.822 


0.829 


Pc (g/cm^) 


86.5 


86.4 


86.3 


86.2 


86.0 


85.8 


85.6 




5.16 10-3 


5.18 10-' 


5.21 10-' 


5.268 10-' 


5.32 10-' 


5.41 10-' 


5.53 10-' 


A 


69.2 


69.4 


69.0 


68.7 


66.9 


65.5 


63.7 



Table 4. Parameters of our purely radiative star for various rotation rates. 




Fig. 7. Same as in Fig.|6]but with Q. - 0.07. Dotted lines are for 
negative contours, in which there is convective instability. The 
thick line represents the points where A^^ = 0. 
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Fig. 8. Brunt-Vaisala frequency profile at pole and equator when 
Q = 0.07. 
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Fig. 9. Differential rotation. The thick line corresponds to Q = 
0.07, which is the background rotation rate. Faster and slower 
rotation are represented by solid black lines and dotted white 
lines, respectively. 



tion in ( fTOl i. together with the equation of state, yields the baro- 
clinic balance, namely, 

^nJVTxV In P\. 

oz 

Now, if we consider the latitudinal variation of the temperature 
along an isobar, replacing the radial variable by the pressure, we 
may write T = T{P, ff), and thus 



dP \ do jp 



We can therefore rewrite the baroclinic balance as 



'dz 



n, ldT\ d\nP 
'dej^ dr 



(34) 



5.2. Differential rotation and circuiation 

We now turn to the dynamical state of our radiative "star". Since 
we simplified the solutions by only focusing on the axisymmet- 
ric steady ones, the dynamics of the star is described by its dif- 
ferential rotation and meridional circulation. We give in Fig. |9] 
the typical shape of isorotation lines in a meridian plane. These 
do not change much with rotation. The main property of this 
baroclinic flow is that the equator rotates faster than the pole. 

This shape of the differential rotation can be understood from 
the torque balance in the momentum equation. Indeed, neglect- 
ing nonlinear and viscous terms, the curl of the momentum equa- 



Since in any direction, pressure decreases with r and isotherms 
are more spherical than isobars (temperature decreases on an iso- 
bar from pole to equator), we find that ^ < 0. The baroclinic 
torque thus makes the polar region rotate slower than the equa- 
torial regions. 

This result is (unfortunate ly!) at odds with the solution of the 
Boussinesq model devised in iRieutordI (l2006h . where the polar 
rotation turned out to be faster than the equatorial one. In the 
Boussinesq model, the baroclinic balance leads to 

— ^ = A^sm^cose, 
oz 
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Fig. 10. Relative differential rotation on the last isobar. 



which shows that rotation increases with z, meaning with lati- 
tude in the northern equatorial region, when the squared Brunt- 
Vaisala frequency is positive (i.e. with stable stratification). 

Hence, compressibility changes the torque drastically. We 
traced this difference to the fact that, in the Boussinesq model, 
density is only a function of temperature; therefore, isochores 
are identical to isotherms, while isochores are strongly influence 
by pressure in the compressible case. 

We also examined the dependence of the differential rota- 
tion on the viscosity and find that it is independent of viscosity 
in the asymptotic regime £ «: 1 . This property is also shared by 
the Boussinesq model, which gives the explanation. As shown in 
iRieutordI (l2006h . viscosity is necessary to determine the differ- 
ential rotation since the inviscid case is degenerate: an arbitrary 
function of r sin (a geostrophic ffow) can always be added to 
u^. Viscosity lifts this degeneracy by "choosing" the geostrophic 
flow such that no tangential stress is exerted on the outer bound- 
ary, while the meridional flow exactly matches the Ekman pump- 
ing. As Ekman pumping and meridional flows are both propor- 
tional to viscosity, this quantity simplifies and we are left with a 
differential rotation independent of it. 

Furthermore, an interesting property is shown in Fig.fTO] ex- 
cept for rotations faster than Q - 0.03 ( ~ 36% of the critical an- 
gular velocity), the relative surface differential rotation keeps the 
same profile. Poles are 10% slower than the equator If the shape 
of this profile is reachable through observations, a comparison 
to this shape will give interesting information on the dynamical 
state of the fluid (presence of turbulence, magnetic fields, etc.). 

We also note from Fig.[TO]that at a high rotation rate, a kind 
of equatorial jet develops between latitudes ±20°. If confirmed 
by more realistic models, such a jet may play an important role 
in mass loss for stars rotating close to the breakup limit. 

The case of meridional circulation, displayed in Fig. (TT] 
shows three cells of circulation. This feature of the meridional 
flow is controlled by the distribution of the Brunt-Vaisala fre- 
quency and the Prandtl number as shown in the Boussinesq case. 
As far as the turnover time associate d with thi s circulation is con- 
cerned, we refer to the discussion in IRieutordI 12006 ) : indeed, ex- 
cept the profile of density, there is no significant change, because 
the turnover time remains controlled by viscosity (probably tur- 
bulent) and the Prandtl number. 

Just as in the Boussi nesq case questi ons of stability would 
require a full study. As in IRieutordI (l2006l) . we nevertheless take 




Fig. 11. Meridional circulation. Solid lines represent counter- 
clockwise circulation and dotted lines clockwise circulation. 
Ps = 10-^ E = 10-\ Q = 0.07, a = 15.2, ctt = 1.5. 




Fig. 12. Directional derivative of along equal entropy lines, 
i.e. (n ■ V)Lj; where n is a unit vector tangential to the isentropy 
surface and directed towards the equator 



a look at the case of axisymmetric barotropic stability. We thus 
compute the angular momentum distribution along an isentropic 
surface; more precisely, we verify that, along such a surface, a 
fluid parcel has an increasing angular momentum when its posi- 
tion gets farther from the rotation axis. This is a necessary con- 
dition for axisymmetric stability, and Fig. [12] shows that this is 
indeed the case. 

5.3. About the effects of the container 

Finally, one may wonder what the role of the box containing 
our "star" is. Obviously, as far as the hydrostatics is concerned, 
the main effect of this box is to reduce the radius (the container 
cuts out the envelope). One may also wonder about the conse- 
quences of the conflicting spherical geometry of the box and 
the spheroidal shape of the mass distribution, on the latitudinal 
flux variations. We find that our model gives a milder ratio of 
the polar to equatorial flux than the von Zeipel model; actually, 
this difference may be reduced when the container is eliminated 
because its spherical geometry may soften the latitudinal varia- 
tions. 
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We also conjecture that the effects of the container on dy- 
namics are fairly weak. Indeed, we note that, while changing Q 
and thus the shape of isobars and isotherms, the form of the dif- 
ferential rotation remains very similar. As the flow comes essen- 
tially from the mismatch of isobars and isotherms, which both 
change with respect to the sphere as Q varies, its invariance at 
low Q pleads for a robust result, independent of the walls of the 
container. 

6. Conclusions 

In this paper we presented a physically self-consistent model of 
a rotating self-gravitating perfect gas heated by nuclear reactions 
and enclosed in a container with a constant absorption coefficient 
radiating like a black body; heat transport in the gas is purely 
radiative, insured by a power law opacity. 

This is a simplified model for a completely radiative star that 
rotates at a significant fraction of th e breakup angular velocity. 
We thus generalized the approach of lRieutordI (2006) by taking 
more realistic properties of a stellar plasma into account. Such 
a fluid is never in hydrostatic equilibrium and a baroclinic flow 
develops. Taking viscosity into account, we computed such a 
flow in its steady state, namely the differential rotation and the 
meridional circulation. Interestingly enough, the results show 
that equatorial regions are rotating faster than polar ones and, 
provided that rotation is less than ~36% of the breakup velocity, 
the surface profile 6Q./Q. is independent of both viscosity and 
rotation rate. 

Although our model is still preliminary, we had a look at the 
anisotropy of the emitted flux to find that the pole/equator ratio is 
less contrasted than what is predicted by the von Zeipel model. 
We also showed that, when rotation is fast enough, equatorial 
regions become convectively unstable. 

In parallel to constructing this model, we tested - for the first 
time to our knowledge - the use of spectral methods in the stellar 
structure computations. Our first results are promising in terms 
of precision and robustness, since we recover the precision of a 
finite difference type model with ten times less grid points. Such 
a gain in efficiency is necessary for all future computations of 
the evolution of a rotating star with a two-dimensional model. 

The next step is of course to confirm these results using a 
model with coordinates adapted to the spheroidal shape of the 
isobars. In this case boundary conditions can be cleanly applied, 
and the surface pressure can be decreased to very low values 
adapted to matching an atmosphere model. Of course, includ- 
ing modeled convection zones is also a necessary step for the 
completeness of models to be used in the interpretation of obser- 
vations. 
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